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INTRODUCTION 


The objective of this study is to develop a two- and three- 
dimensional comprehensive spray combustion computer simulation 
code to study the complex physical processes involved in a 
liquid-fueled combustor. Efforts of this study in the last six 
months include: ; 1) implementing the Stochast ic/two-equat ion 

model for turbulent droplet dispersion calculation; (2) testing 
the transient particle tracking methodology in a benchmark flow 
field; (3) upgrading and test the two-equation k-E model for the 
ARICC code. In the following, progresses in these studies will 

be briefly described. 

(1) implementations of the Stochastic Lagrangian particle 
tracking model into the MAST code have been completed. We 
utilized the operator-splitting technique such that the two-way 
coupling between the two phases is accounted for in a multi- 
corrector procedure. This operator-splitting technique eliminate 
the global iteration processes used in the conventional 
SIMPL/PSIC (particle source in cell) method. This method is time- 
accurate, and has been shown to very efficient for transient 
spray calculations. 

(2) Testing of the above mentioned two-phase methodology 
for single injector spray has been carried out. The testing 
conditions for the Horiyasu's experiment setups are shown in 
Table 1. At the injection exit, the distribution of droplet 
sizes is modeled by a X-square function. Particle/droplet 
interaction was modeled by random sampling from assumed 
probability distribution of flow turbulence. For the k-E model 
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used, this probability distribution is Gaussian. Figure 1 and 2 
show the penetration of spray at different chamber pressures. 
Comparisons of the computed spray penetrations with experimental 
data are shown in Figure 3. The agreement is excellent and shows 
the time-accuracy capability of the current methodolgy. The 
efficiency assessment of the MAST calculations is summarized in 
Table 3. It can be seen that the MAST reduces CPU time by one 
order-of-magnitude and very much particle-number insensitive. 
The robustness of this methodology is demonstrated for this test 
case. Further testings are underway for hollow-cone sprays. 

(3) Turbulence model upgrade for ARICC and the testing of 
the ARICC performance on benchmark turbulent flows have been 
completed. Due to the explicit ALE-ICE method used in the ARICC 
code, the inclusion of the k-E model highly reduced the time-step 
for marching solutions. The efficiency of the ARICC code with 
advanced turbulence models was seen to deteriorate significantly. 
Detailed implementation and validation studies can be found in 
Lee's Master Thesis which is funded by this grant and is enclosed 

as an appendix. 
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TABLE 1 


SINGLE- ORIFICE INJECTION PARAMETERS(Horiyasu) 


Chamber 
Gas Pressure 
{ATM) 

Injection 

Velocity 

{m/ sec) 

Gas 

Density 

(kg/m 3 ) 

Mass 

Flow 

(kg/ sec) 

Sauter Mean Eddy 
Radius(SMR) Viscosity 

fjm (m 2 / sec) 

1 

122.2 

1.123 

0.00726 

5.0 

7.1 x 10“ 

30 

102.5 

33.70 

0.00609 

5.0 

5.0 x 10“ 


Fuel: Diesel fuel oil, p s = 840 kg/m' 


Ambient Gas: Nitrogen 
Nozzle Diameter: 0.3 mm 



TABLE 3 


EFFICIENCY ASSESSMENT (CPU Time) 



SINGLE-ORIFICE SPAR 
41 x 61 Grid 
300 Time Steps 

HOLLOW CONE SPRAY 


MAST-2D 
Particle # 

600 126.9 sec 

1200 135.7 sec 

400 74.9 sec 


TEACH/PSIC 
Particle # 

800 1420 sec 


31 x 31 Grids 
200 Time Steps 


800 934 sec 


1000 


88.3 sec 
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ABSTRACT 


The purpose of this study is to validate and further develop an existing Compu- 
tational Fluid Dynamics code for simulating complex turbulent flows inside a liquid 
rocket combustion chamber. The ARICC (Advanced Rocket Injector/Combustor 
Code) Code is simplified and validated against benchmark flow situations for lami- 
nar and turbulent flows. The numerical method used in ARICC Code is re-examined 
for incompressible flow calculations. For turbulent flows, both the subgrid and the 
two equation k - e turbulence models are studied. Cases tested include idealized 
Burger’s equation in complex geometries and boundaries, a laminar pipe flow, a high 
Reynolds number turbulent flow, and a confined coaxial jet with recirculations. The 
accuracy of the algorithm is examined by comparing the numerical results with the 
analytical solutions as well as experimented data with different grid sizes. 
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CHAPTER I INTRODUCTION 


Combustion flows involving liquid propellants inside a liquid rocket engine rep- 
resent one the most complicated engineering flow systems in operation. The com- 
plexity stems from the existence of multiple zones with very different time and 
length scales within the same physical domain (Figure 1.1). The zones are in close 
proximity to each other and the processes are usually strongly coupled. In the 
injection zone, for instance, the major physical characteristics include those of mul- 
tiphase flows, recirculation, high shear stresses, and steep pressure gradient. The 
multiphase flows nvolved are : a multispecies gaseous phase, incompressible liq- 
uid phase, and a ^articulate droplet phase. The dominant processes are liquid jet 
breakup and atomization, droplet-gas interactions and fluid-wall interactions. A 
short distance away from the injector, the characteristics of the flow become those 
of steep concentration gradients around droplets, turbulent mixing and diffusion, 
and dilute gas-liquid suspension flow. The critical processes become those of droplet 
transport, evaporation, and droplet heatup. Further downstream, the fuel/oxidizer 
mixture begins to combust. The flow is now characterized by high heat fluxes, a rel- 
atively large number of gaseous species, and temperature extremes that necessitate 
the consideration of real gas properties. The key processes are kinetic and equilib- 
rium chemical reactions and species diffusion. Finally, in the post-combustion or 
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expansion zone, the fluid properties range from a “frozen” composition to near “dy- 
namic equilibrium". There is high velocity gradient near walls and at the nozzle. 
The dominant processes of interest are transonic/supersonic flows and boundary 

layer effects. 

As shown in Figure 1.1, all zones are coupled with each other. The first three 
zones in a typical rocket engine preburner, for example, the Space Shuttle Mam 
Engine (SSME), comprise a volume of only about 1-inch diameter by 2-mch length. 
While different dominant processes can be identified in different zones, they are 
intrinsically coupl0d to each other. 

In 1985, a comprehensive model for the simulation of detailed three-phase 
combustion flows inside combustion chamber was developed by Rocketdyne. The 
Rocket dyne Code (ARICC) is developed from CONCHAS-SPRAY of LANL (Los 
Alamos National Laboratory). The original CONCHAS-SPRAY code solves the 
equations of transient multicomponent chemically reactive fluid dynamics, together 
with those for the dynamics of an evaporating liquid spray. With few exception, 
reactive flow problems of practical interest are far too complex to be solved ana- 
lytically. Quantitative theoretical analyses therefore require the use of numerical 
methods. CONCHAS-SPRAY is a time-marching code that solves finite difference 
approximations 1o the governing differential equations. The transient solution is 
marched out in a sequence of finite time increments called cycles or time steps. 
Values of the dependent variables on each cycle are calculated from those on the 
previous cycle. CONCHAS-SPRAY is a two-dimensional code, which assumes 
that the dependent variables depend on only two of the three spatial coordinates 
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because of symmetry. The effects of turbulence are represented by a simple subgnd 
scale (SGS) turbulence model, whose use is optional. The option is also provided 
to calculate boundary layer drag and heat transfer from a modified turbulent law 
of the wall. CONC HAS-SPRAY utilizes a partially implicit numerical scheme that 
is a variant of the ICE method (Harlow and Amsden, 1968, 1971). Spatial differ- 
ences are formed with, respect to a generalized finite-difference mesh or grid, which 
subdivide the region of interest into a number of small quadrilateral cells, or zones. 
The corners of the cells are called the vertices. The position of the vertices may be 
arbitrarily specified as function of time, thereby allowing a Lagrangian. Eulerian, 
or mixed description. Since the locations of the vertices are arbitrary, the cells are 
arbitrary quadrilaterals. This type of mesh is called an ALE (arbitrary Lagrangian- 
Eulerian) mesh (Hirt, Amsden and Cook, 1974; Pracht 1975), and is particularly 
useful for representing curved and/or moving boundary surfaces. Evaporating liquid 
sprays are represented by a discrete- particle technique (Dukowicz 1980), in which 
each computational particle represents a number of similar physical particles. 

Simulation of gas phase and particle phase was completed in CONCHAS- 
SPRAY. In ARICC code, the volume of fluid method was added to account for 
dense spray (Figure 1.2). While the liquid phase is assumed to be incompressible 
with a constant ciensity and temperature. The “pseudo-incompressible treatment 
was used in ARICC, the inclusion of limited compressibility made the liquid iter- 
ation process much more stable and allowed much faster convergence (Liang et al. 
1985). 

The ARICC code is designed to simulate two-phase multi-species reacting flow 
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with spray droplets in either a two-dimensional or an axisymmetric with swirl con- 
figuration. The code is real-time accurate. It accounts for all the important physical 
processes involved in commonly encountered combustion devices either vigorously 
or using empirical models. Radiation, which is not significant in the case of liquid 
rocket combustion chambers, is not modeled (ARICC User s Manual). 

The Deardorfi’s subgrid model is an algebraic model in which the “mixing 
length” is prescribed by the local grid size. In the last two decades, the advances 
of turbulence modeling has reach the point where reasonable good predictions of at 
least the mean velocity fields can be made by some advanced models. The purpose 
of this study is incorporate the so-called two equation kinetic energy (k)- energy 
dissipation rate («) turbulence model into the ARICC code for better presentation 
of turbulence field simulation. The k - e model accounts for the history effect, 
as well as diffusion, sources of the characteristic velocity and length thus is more 
capable of handling complex flows such as the ones encountered in typical liquid 
engines. In the following chapters, numerical aspects of ALE-ICE schemes will be 
described. Several benchmark testing cases for laminar as well as turbulent flows will 
be performed to validate the implementations and performances of the turbulence 
models and the ALE-ICE scheme. Finally, conclusion and recommendations will be 


made. 



IStON 


COMBUSTION 

zone 


EXPANSION 

ZONE 


Coaxial Injector /Combustor 




V O F (Volume of Fluid) 


• ACTIVE CELLS DISTINGUISHED BY THEIR FLAG VALUES: 



FULLY LIQUID PARTIALLY LIQUID FULLY GAS 

F = 0.0 F = 0.7 F = 1.0 


NET GASEOUS VOLUME (EXCL. DROPLETS) 

DEF F ' LIQUID VOL + NET GAS VOL. 

• DROPLETS MAY EXIST IN CELLS WITH F > 0.0 

• SINGLE-VALUED THERMODYNAMIC PROPERTIES FOR ALL CELLS 

• PRESSURE IN FULLY LIQUID CELLS DETERMINED BY PRESSURE 
ITERATION SCHEME FOR IMCOMPRESSIBLE FLUID 


Figure 1. 2 Volume of Fluid Algorithm 




CHAPTER II MATHEMATICAL MODELS 


A. Cnvprning Equations 

The governing equations in axisymmetric fluid flows are shown in cylindrical 
coordinate with x-y-9 coordinates and u-v-w velocities. The x direction is the axial 
coordinates and the y direction is the radial coordinates. 

For conservation of mass, the governing equation is 

% + 1 y& pyu)+ T y {py ' l)i=0 (21 ' a) 

For conservation of linear momentum of Newtonian fluid, the governing equa- 


tions are given by 


l (pu) + i( ^ = -|H + i[|.(»T„) + %-(yr„)) (2.2 .a) 

dt KP ’ y dx dy dx y l dx dy 


d 1 .dypuv dypvv dp l,d d , ^ 

+ -aT J = “a , + i l fe ( » T " ) + dy {yT ” >] 


V 


Tee | pw 

y y 


(2.2 .b) 


The angular momentum equation, which determines the swirl velocity w, is 
given by 

j t (ypw) + i ^ y2 pwu) + §^ {y2 pwv)] = (2 ' 3) 

For conservation of internal energy, the equation is given by 

d d d 

■SpD + + si iypIv)] = ~^ l di {yu) + T y {yv)]Jr 


df 


y l dx 
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where 
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J X m dx 3 J 
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for the energy equation, the heat flux is 


du v Ov 

Ox y Oy 


OT 

J ' = - K W 

r .,dT. 

J y = ~ K % ) 

where p is the density, p is the pressure, u is the velocity in axial direction, t- is 
the velocity in n dial direction. T is the temperature. I is the internal energy of the 
fluid flow, p is the viscosity of the fluid and K is the thermal conductivity. 

For hot gases, equation of state is described by the ideal gas law 


p - pRT 


(2.5) 
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where R is the un versal constant. 

The preceding equations have heen ffi'^n in forms appropriate for laminar flows. 
In order to treat turbulent flows, the equations must be suitably averaged. Ensemble 
averages are used here. Whatever type of averaging process selected, instantaneous 
dependent variab.es are separated into mean and fluctuating components. The 
averaged equatior s are considerably simpler if the mean values are mass weighted 
(the so-called Favre averaging procedure) for compressible turbulent flows. 

In order to slow the turbulent motions and their interactions with the mean 
motion more exp icitly, we follow a procedure due to Osborne Reynolds. Let us 
apply this procecure to the flow of a. incompressible fluid with constant \iscosit\, 
constant density, isothermal flow. In this approach the instantaneous quantities are 
decomposed into mean and fluctuating parts, i.e.. 


u j = l t + u 7 
p = P + P 

After insertion ir to equations (2.1) and (2.2) the following time averaged equations 


result 


■«± + A = ” + ^ 

of 3 dxj d.r , dxj dxj dx t 


16) 


For an incompressible fluid and constant viscosity, we come to the conclusion that 
the turbulence terms pu'u' can be interpreted as stresses on an element of the fluid 
in addition to the stresses determined by the pressure P and the viscous stresses. 


Because Reynolds was the first to give the equation for turbulent flow, the turbulent 


stresses pu\u^ are often called Reynolds stresses. 
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In comparing the turbulence stresses in the equations of motion with the cor- 
responding stresses caused bv viscositv effects it is tempting to assume that the 
turbulence stresses act like the viscous stresses, that is. that they are directly pro- 
portional to the velocity gradient. Thus 


dUi dUj . 

+ ii‘ ) 


( 2 . 7 ) 


This assumption was made by Boussinesq, who introduced the concept of an appar- 


ent, or “turbulence,'’ or “eddy” viscosity /q. According to Boussinesq's concept, 
the eddy viscosit} /q has a scalar value. 

In this approximation the averaged flow equations become identical in form to 
the laminar ones; the transport coefficients (i.e., viscosity and thermal conductivity ) 


are simply replaced by the appropriate turbulent transport coefficients which are 
much larger than the laminar values because of the additional transport caused by 
the turbulent fluctuations. Therefore, we use the equations summarized above even 
when the flow is turbulent, but with turbulent contributions added to the laminar 
values of the transport coefficients. The transport coefficients in ARICC are thus 


taken to be 


/< = HL + 


( 2 . 8 ) 


K = HC v /P r 

In this study, the dependence of m on temperature is modeled by a curve fit, fJ-L 
is laminar viscosity, P r is the Prandtl number based on c„, c„ is specific heat at 


constant volume, fit is a turbulent viscosity. 
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The eddy viscosity fi t thus defined becomes a property of the flow encountered 
?ud requires modeling. The hierarchy of turbulence models has been reviewed 
extensively in the literature (see Deardorff (1971) and Jones and Lauder ( 1972)). In 
this study, the original zero equation subgrid turbulence model originally used m the 
CONCHAS-SPRA Y and ARICC code will be validated against incompressible flows. 
In addition, a two- equation turbulence model ( k - e model) will be incorporated 
into the ARICC code. Brief descriptions of these two models are given below. 


B. The Subgrid Turb ulence Model 

The subgrid model was suggested by Deardorff (1971). The local Reynolds 
stresses which arise from the averaging process were simulated by an eddy coefficient 
with magnitude limited in some way by the size of the averaging domain. When 
this domain is considered to be the grid volume of a detailed numerical integration, 
the eddy coefficient „ ( becomes a “subgrid scale" of “SGS" eddy coefficient. The 
formulation, which allows n t to be variable in space and time, is formulated as: 

= ±pK 2 D A 2 (D 2 u + D 2 22 + DU + 2 D\ 2 + 2Dl , + 2T> 2 2 3 )’ (2-10) 

v2 


where 


_ du n dv 

Du=2 - D 22 = 2-, D„ 


9 — 


du dv u / w \ n -,,^-fi) 

Di 2 = 7 T + 7r: D 13 = y «-(-)> D 23 y q ( ) 


d,w 


d , w 


dy dx 


dx y 


dy y 


( 2 . 11 ) 


where A is a repi esentative grid interval and was taken to be the largest side length 
of the regular ceil. D t] is the magnitude of the local velocity deformation calculated 
on the finite-difference grid, and K D is the dimensionless constant. Lilly (1967) has 
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estimated I\ D wrh approximate knowledge of the Kolmogorov inertial-subrange 
constant a He assumed that isotropic turbulence with such an inertial subrange 
was present in the problem being simulated numerically on scales both much greater 
and much less then the grid interval A. He also assumed Ax = Ay = A;: — A, and 
took the Reynold ; averaging volume to be A 3 . He found that 

Kq = 0.23a - ^ C 2 - 12 ) 


The essential tec! niques recommended by Williams (1969) were used independently 
by Deardorff (1970) in a steady of turbulent channel flow at large Reynolds number. 
I\ D was suggested to be equal to 0.17. 


C. The k - e Turbulence Model 

The field of turbulence modeling for single-phase flows is a rapidly expanding 
one and many p oposals have been suggested. The chosen two-equation model is 
used to generate a turbulence length and velocity scale and these will be used to 
form a (non-constant) eddy viscosity. An appropriate velocity scale is u' = fc* • The 
k — e model adopted here is the one developed by Jones and Launder ( 1972). In 
their paper, the equations was given by 


Cpk 1 , d 


1 f d u dk 


dt 


+ + g-(vM-)] = 


y ' dx XJr " ~ ’ ' dy 

d /i dk 


n ( — V^")} + G-c D pe 

dy ok ay 


dpt 1 . 
dt ^ y dx 
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1 r d p dt 


+ -\—a<JP u ^+Ty [lll,Ve]] = y { d~x ( V,' J di' 


d y dt (c\ e tG — C 2 t pt~ ) 

+ ^W + 5 


(2-13) 


(2.14) 
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where 


fin 


fiv 


fin fi 


fi 


fi 


^ {21(77-) +( o-r + <->' 


4- \ 


'I' + w 


fi r ' d” ‘ ;i ' ■ 'v 

is then expressed in terms of the turbulent kinetic energy k and the turbulent 
energy dissipation rate e via the relation 

(2.15) 


fi 


t — Cupk /e 


Philosophically, the strongest motivation for turning to more complex models 
is the observation that the algebraic model evaluates the turbulent viscosity only in 
terms of local flow parameters, yet a turbulence model ought to provide a mechanism 
by which upstream effects can influence the turbulence structure (and viscosity) 
downstream. Further, with the simplest models, ad hoc additions and corrections 
are frequently reemred to handle specific effects, and constants need to be changed 
to handle different classes of shear flows. To many investigators, it is appealing to 
develop a model general enough that specific modifications to the constants are not 

require to treat different classes of flows. 

Note that the equations (2.13 ~ 2.15) used above are valid only for high 
Reynolds number turbulent flows. In order to simulate low Reynolds number tur- 
bulent or transition flows including the near wall flows, the high Reynolds number 
k - e model must be modified to incorporate viscous and low Reynolds number 
effect. In this thesis, only high Reynolds number k — e model is considered. For 
this reason, the wall boundary conditions have to be constructed to avoid sublayer 
viscous effects and the wall function approach (to be described in Chapter 4) is 


adopted in this study. 
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In the equation (2.13) ~ (2.15) the quantities <r k , *<, cu, c 2t , and c D are 
mo^el constants. Following Jones and Launder ( 1 Q72) *be wflues used are hstmg in 

Table 2.1. 


Table 2.1 Constant values for k - e model 


db cr f c.\ e c<it cd 

1.00 1.30 1.44 1.92 0.09 1.00 


Numerous other two-equation models have been suggested, Rubesin (1977) 
shows several comparisons between these models for incompressible flow, and overall 
they perform quite well, it is difficult to identify the best from the comparisons he 

has shown. 

The k-e turbulence model has enjoyed wide use because of its ability to predict 
the mean velocity field and spreading rate of many turbulent shear flows. A recent 
review of the applications of this model to a wide range of problems is given m Rodi 
(1982). Application of the k - e model to heat and mass transfer problems is given 

in Shih (1982). 

Despite the enthusiasm which is noted from time to time over two-equation 
models, it is perhaps appropriate to point out again the two major restrictions on 
this type of model. First, two-equation models of the type Boussinesq approxi- 
mation holds. Iii algebraic models, fi t is a local function whereas in two-equation 
models fi t is a more general and complex function governed by two additional PDE’s. 
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If the Boussinesq approximation fails, then two-equation models would fail. Obvi- 
ously. in manv flow- the Boussinesq approximation models reality rloselv enough 

for engineering purposes. 

The second shortcoming of two-equation models is the need to make assump- 
tions in evaluating the various terms in the model transport equations, especially in 
evaluating the third-order turbulent correlations. This same shortcoming, however, 
plagues all higher order closure attempts. These model equations contain no magic, 
they only reflect he best understanding and intuition of the originators. We can 
be optimistic, however, that the models can be improved by improved modeling of 

these terms. 



CHAPTER III NUMERICAL METHODS - 

THE ALE-ICE SCHEME 

The numerical method used in the ARICC code is the ALE-ICE (Arbitrary 
Lagrangian-Eulerian Mesh, Implicit Continuous-fluid Eulerian Technique) scheme 
which utilizes the fractional time step concept and solves the governing flow equa- 
tions in unsteady forms. The temporal domain was discretized into time steps 

At = t n+i -t n ( i=0,l, 2,...). 

From t n to was called a cycle. There are two phases to be performed in 

one cycle. From t' 1 to the first intermediate step was called phase A. From phase 
A to next intermediate step was called phase B. Then from phase B. we go to next 
step phase C which is equivalent to f' i+1 step. Phase A is an explicit Lagrangian 
calculation, phaie B is an implicit pressure correction and phase C or n+1 step is 
a rezone calculation. 

This separation of a calculational cycle into a Lagrangian phase and a con- 
vective flux, or rezone, phase originated in the Particle-in-Cell numerical method 
(Harlow 1955, A nsden 1966). and has since been used in many hydrodynamic com- 
puter codes. Ir the present technique the different phases can be combined in 
various ways to suit the requirements of individual problems. For example, in high 
speed problems, in which the Courant stability condition is not likely to be violated, 
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an explicit calculation is acceptable and the phase two iteration may be omitted, 
and for an explicit Lagrangian calculation only phase one used. 

The scheme uses arbitrary Lagrangian-Eulerian Mesh. This type mesh is made 
up of arbitrary quadrilateral. The equation was discretized by control volume 
or integral-balance approach, which preserves the local conservative properties of 
the differential equations. Spatial calculation domain is divided into a number of 
nonoverlapping control volumes such that there is one control volume surrounding 
each grid point. The differential equation is integrated over each control volume. 
The grid vertices may move in an arbitrarily prescribed manner. This capability 
includes the Lag] angian and Eulerian description as special cast s. 

A regular cell is shown in Figure 3.1. The index of the regular cell (ij) were 
regarded as horizontal and vertical coordinate. The indices (ij) also label the ver- 
tices, with the understanding that vertex (ij) is the (logical) lower left corner of cell 
(ij). The "center ( rjj . ) of the cell is defined by 

x c n = 7(^1 + x 2 -|- i 3 + x 4 ) 

Vri = 7 ( 1/1 + V2 + V 3 + Vi ) ( 3,:l ) 

J 4 

In general, the point (x^y^) is not the center of mass or volume of cell (ij). 

The area or’ the cell in Figure 3.1 following CONCHAS-SPRAY (1982) and 
KIVA-II (1989) is calculated by the algebraic formula, 
area = trial lgle 123 + triangle 134 

\TR= i(yi(x 2 -x 3 ) + y 2 (x 3 -xi) + y 3 (xi -x 2 )) 
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ABL = \(yi(x 3 - x A ) + y z {x 4 — an) H- y 4 (*i - * 3 )) 

A ti = ATR + ABL (3-2) 

The volumes < orresponding to triangles 123 and 134 may be shown to be 


/ TR = f ydxdy = — (j/i + J /2 + ?73 )ATR 
J A( 123 ) ^ 

VBL = f ydxdy = -(yi + y 3 + y 4 )ABL (3.3) 

J A(13 4 ) 3 

The total cell volume is then given by 


Vij = / ydxdy — VTR + VBL (3.4) 

J-Uj 

Momentum cell (ij) is centered about vertex (ij), as shown in Figure 3.2. A 
momentum cell hr s four of its corners at the centers of the four associated regular 
cells, the other four at the midpoints of the regular cell sides which meet at vertex 
(ij). Momentum 'ell's area and volume are M”) and V ” 1 . Those calculations are 
similar to the regular cell’s calculation (equation 3.1 ~ 3.4). 

If the finite-difference representation for the equations has the conservative 
property, we musi establish that the discretized version of the divergence theorem 
is satisfied. The cifferential equation is integrated over the area of a typical cell or 
momentum cell ai d use divergence theorem to transfer area integral to face integral. 


VF dA 


A 


-l 


F • nds 


We normally check this for a control volume consisting of the entire problem do- 
main. To do this the integral on the left is evaluated by summing the difference 
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representation of the equation at all grid points. If the difference scheme has the 
conservative property, all terms will cancel except those which represent fluxes at 
the boundaries. It should be possible to rearrange the remaining terms to obtain 
identically a finite- difference representation of the integral on the right. The result 
will be a verification that the mass flux into the control volume equals the mass 
flux out. When the governing equation can be written in divergence form, we can 
be guided in this process by employing the Gauss divergence theorem to obtain 
the correct mathematical formulation for the physical law for a control volume. In 
practice, the control volume method has a history of leading quickly to expressions 
that prove to be more accurate than other possibilities near boundaries, probably 
because the method keeps the discrete nature of the solution method in view at 
all times. The distinctive characteristic of the control volume approach is that a 
“balance” of some physical quantity is made on the region in the neighborhood of a 
grid point. The discrete nature of the problem domain is always taken into account 
in the control volume approach which ensures that the physical law is satisfied over 
a finite region rat ier than only at a point as the mesh is shrunk to zero. It would 
appear that difference equations developed by the control volume approach would 
almost certainly have the conservative property. 

Spatial differences are performed by integrating the differential term in question 
over the area of a typical cell (or momentum cell). If the term is of gradient 
or divergence form, its area integral is usually converted into a surface integral 
using the divergence theorem. The area integral of a time derivative is related to 
the derivative of the integral by means of the Reynolds transport theorem. The 
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area and surface integrals are performed under the assumption that the dependent 
variables with control volume approach are assumed uniform in each cell and/or on 
each face. Velocities are set at vertex and are regarded as uniform in momentum 
cell. All quantities are regarded as uniform with the value at the center of the cell 
face. Then, the line integrals over cell face are approximated by sums over the cell 
faces as shown in Figure 3.2. 

/ F • nds = Fp, • n a A s Q (3.5) 

Js Q 

where F 0 is the valv e of F at the center of cell face a, n a is the outward unit normal 
to the cell face, As £l is its length, and the faces are numbered in a counterclockwise 
order. 

n a As a — 1q X k = /aryl ZarxJ (3-6) 

where 1 Q is the ve:tor of length along cell face a in the counterclockwise 

direction, k is the unit vector out of the plane, bidav are the y components of 
\ a (Figure 3.3). 

The velocities are located at vertex, then 

Uj j — u(xjj,yij) 

Vij = v{x t j,y tJ ) ( 3 ' 7 ) 

where u,v are the x, y components of fluid velocities. Thermodynamics variables 
and other scalar variables such as turbulence kinetic energy and energy dissipation 
rate are located at cell centers. 


Qxj = Q(x c ij,y c ij) 


(3.8) 
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where Q = p, p, T , I, fc, e. 

A. PViasp A : The Explicit Lagrang ian Calculation 

The explicit Lagrangian calculation was used in phase A. In this step velocities 
are advanced explicitly in time using pressure gradients and body forces computed 
from the currentlj available pressure and mesh coordinates. If viscous, elastic or 
other stresses are desired, they may be included at this stage as well. The total 
energy of each cell is also advanced in time to account for the work done by the 
body forces and other stresses, except those of pressure. Pressure work terms are 
included only after the implicit pressure calculation in phase B. This delay permits 
time-advanced pressures to be used in computing the work and ensures consistency 
with the velocities coming out of phase B. The vertices were assumed to move 
locally with the velocities. From this assumption, the convection terms of transport 
equations vanished. From continuity equation, we have 

= f ypdA = p'ljVt] (3-9) 

JAij 

From the Lagrangian assumption, the vertex position can be calculated by 

xt, = X*ij + < A t 

y?i = y." + Ai (3 - 10) 

and the volume I'd can be calculated by the equations (3.1) ~ (3.4) with all x”’s 

l J 

and y n, s replacec by x A 's and y A 's. The p A j, p A were calculated by 

pti = M*!V* 


(3.11) 
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(3.12) 


The momentum equation was discretized by integrating over the momentum 
cell area A*. In momentum cell the velocity is regarded as uniform, but the 
density was calculated by the associated four regular cell’s density. 


Therefore, 


[ ypuijdA = uij[ ypdA — MijUi 

Jay; 

[ ypVijdA = Vijf ypdA = MijVi 

Jay; Ja ?, 


(3.13) 


where M>, is the mass of momentum cell (ij). The momentum equation in phase 
A was discretized by 


n l J 


Z = _y» pi l n ax + ^(yr y „ Ax + yr yx Ay - T M A m )Z + M% ^ 


a= 1 


a=l 


MMsddMl = £ p ; i; f + £(</n„,Ax + »r„Ay); (3.14) 

At 0*1 0=1 

where the p n a is the pressure of the regular cell in which face a lies, the summation 
are over the four faces which are orthogonal to the pressure directions. Ax and Ay 
are the distance of the cell center to the cell face in x direction or y direction. 

The diffusion terms summation are over the eight faces of momentum cell (ij), 
it has been show i in Figure 3.2. For evaluating viscous stresses, the regular cell 
was partitioned into four quadrants by the momentum cell faces which meet at the 
cell center 0 (Figure 3.4). These faces are labeled R, T, B, L. Each vertices has two 
faces from one regular cell, and eight faces from four associated regular cells. Faces 
B, R are two of he eight faces of momentum cell centered at vertex 1, faces R, T 
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are two of the eight faces of the momentum cell centered at vertex 2 and so on. The 
velocities at point 0 are defined as the average velocities of the four vertices values. 

Uq = -(«1 + «2 + u 3 + u i ) 

4 

V 0 = ^(uj + v 2 + V 3 + t? 4 ) (3.15) 

The viscous stresses on faces R, T, L, B were described below. First, the partial 
derivatives of u, v with respect to x, y in faces R which shown in Figure 3.4 were 
calculated by the triangle rule. 

r du , 2/20^10-1/10^20 

R — " 

l dx Xio 2/20 — x 20 2/10 

.dv _ Zip u 2 q - ^20 ^10 (3.16) 

9y R xio 2/20 - x 20 2/io 

where 

Uab = U a - Ub, V a b = v a ~ v b 
X a b ~~ X a Xbi Uab Va Vb 
a, 6 = 0,1, 2, 3, 4 

and the others faces T, L, B are obtained by the permutations ( R, 1, 2 ) - ( T, 3, 
4 ) _► ( B, 4, 1 ). The diffusion term is evaluated as the sum of four contributions, 
one from each of the regular cell shown in Figure 3.4 contributes to each of the 

vertices 1, 2, 3, 4. 

B. Phase B : Pressu re Correction 

In phase B, the vertex positions are determined by the equations 

+ «g A< 
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vS = v5 + "S A* 


> 1 } 


(3.17) 


From those veitex positions *g,yg, the cell volume V,f in phase B can be 
calculated by the ecuation (3.1) ~ (3.4). The cell mass is not changed in this step, 
because the continuity equation. 


M b = M a 


(3.18) 


Then the density 


P?, = = M{}/V* = P$V4/V 


rB 

>3 


(3.19) 


Therefore, the pressure can be calculated by 


B n rnB 

Pij ~ Pij *ij 


(3.20) 


For isothermal process, we have T® — Tf- — T,". 

In phase B, the momentum equations was reduced to 


uf - ufj 

Vfd-i! lJ 

A t 


^ =-».”E(pf -P'oKy 


a= 1 


M 


a v A z V A 

') At 




a—1 


(3.21) 


In equation (3.21) the summations are over faces of momentum cell (ij), p a is 
the pressure of the regular cell in which face a lies. 

The methodology described above is valid for compressible flows since equation 
of state is utilized. As for the unsteady Navier-Stokes equations, the numerical 
difficulty lies in the constraint V • V = 0. In the steady case, this difficulty can 
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be surmounted by using the so-called Pseudocompressible method. The method 
has been introduced by Chorin (1967). The principle of the method is to consider 
the solution of the steady equations (2.1) as the limit when t - oo of the solution 
of unsteady equations obtained by associating the unsteady momentum equation 
(2.2) with a perturbed divergence equation in order to get a system of equations of 
evolution which can be easily solved by standard methods. 

The Chorin method is established by first writing a perturbed continuity equa- 
tion 

+ c 2 P V • V = 0 (3-22) 

dt 

where c 2 is an arbitrary constant. This equation has no physical meaning before the 
steady state d/dt = 0 is reached. So the constraint V • V = 0 is satisfied at conver- 
gence only. The method, which consists of solving equation (2.2) and (3.22) can be 
called a pseudo-ur. steady method because the time t involved has no physical mean- 
ing. The parameter c 2 in equation (3.22) must be chosen to ensure convergence, 
i.e., to ensure the existence of a steady numerical solution of the system. 

The term “pseudocompressible method” (Chorin 1967) can be derived from the 
Navier- Stokes equations for a compressible fluid whose state law would be 

P = c 2 p 

with c 2 = constant. However, possible numerical difficulties can be associated with 
the use of a very large value of c 2 , and hence the pseudocompressible method will 
likely have the most value in the computation of steady solutions. In this way, it 
can be considered a procedure to build a special iterative method for solving the 
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steady problem. For incompressible flows, the momentum equations must coupled 
with the continuity equation. Continuity equation was discretized by 

(p B - P") + PC 2 (^ - 1) = 0 (3.23) 

The equations (3.17) ~ (3.21) for compressible flow or (3.17), (3.21), (3.23) 
for incompressible flow are an implicit system of equations for unknown quantities 
t b i,b u b and p B One of the major problems in machine computations is 

’ ij y l 3 * 

to find an effective method of solving a system of n simultaneous linear equations, 
particularly if n is large. There is, of course, no best method for all problem because 
the goodness of a method depends to some extent upon the particular system to be 

solved. 

This implicit system can be solved by a Successive Over Relaxation (SOR) 
scheme. We now oroceed to describe the SOR. Let the iteration index be v, which 
will be displayed as a superscript in parentheses. SOR is specified by giving the 
prescription for advancing all quantities from iteration v to iteration u + 1. The 
vertex positions are advanced first. 


x 


<r +1) = x", +u(;> At 


V 


A( 


(3.24) 


the quantities (x<T +1 \ itjf >») then determine V^ + " in the usual way. The pressure 
change is computed by the relaxation formula for compressible flow 


(iH-l) _ „(") fl. .[„("> _ U— M A 1 

Pij 'Pi] ~ PijlPij 

V ij 


(3.25 .a) 
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where 


ry(v) rjiA 


for incompressible flow, the equation changed to 


i „(*0 

rlPij ^yCi'+l) ynj 


</+l) ( v) 


p!;’ - (Jo [<?!;’ 


. 2 / O' 


i)] 


(3.25.6) 


and fa is a rela> ation coefficient defined below. Finally, the velocity is updated 
according to 


.wH* - *" = M W> ~ At »« - p'X» 

a 

Al^!; +1) = calMfjvfj - Aiy” 

a 

Since equation (3.25) determines the pressure change 




It is convenient to rewrite equation (3.26) in terms of 


= -sr 11 




We then obtain 


Z n 

ay 


(3.26) 


(3.27) 


(3.28) 


ifftu If = -A tyf^Sp^r, 

a 

MljSvtf = -A(y; £ (3 ' 29 > 

a 

Equation (3.29) is more convenient for numerical purposes than equation (3.26), 
as it allows one to avoid saving the ufyvfy 
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The iteration is initialized by setting 


,.(o) 4 (o) _ A J 0) = p n . + Ap 

u ,j ~ U i J’ V 'j ~~ o’ "'3 F 


(3.30) 


where Ap is an estimated change in overall pressure level (defined below), and the 
final converged values are given by 


(") „.(*') „") 




,P,n u 


(") ..(•') 






(3.31) 


of course, the limit v -» oo is achieved for all practical purposes at some finite 
value of v\ this value is determine by a convergence criterion and the iteration is 
terminated at that point. After the iteration has converged. V? is known and 
the densities pf} are calculated by (3.17). All quantities are updated by its relative 
equations. 

The relaxation coefficient ft x j for compressible flow is given by 

/% = »[! (3 32.«) 

for incompressible flow is changed to 


, .dSij i 

^ = w{ d^ ] 


(3.32.6) 


where 

yiv+l) 

= CpSJ* - p?>) + — 11 

where V l} and T, t are here to be considered as function of p XJ ; these functions are 
implicitly determines by equation (3.1), (3.2), (3.23), (3.24), (3.29). The derivative 
in equation (3.31) is evaluated numerically by inserting a small pressure perturba- 
tion into cell (i,j) and calculating the resulting change in ( Tij/Vij ) for compressible 
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flow or Sij for incompressible flow. The coefficient w in equation (3.31) is the usual 
overrelaxation parameter, which may be varied in the range 0 < w < 2 to accelerate 
convergence. Best lesults are usually obtained with w > 1 . 

C. Phase C or n+1 Step : Rezone Calc ulation 

In phases A and B, the vertices were assumed to move with the fluid, based on 
the Lagrangian calculation. If the calculation in phase C used Lagrangian technique, 
the rezone procedure is not needed. If the Eulerian technique was used, the new 
mesh should be remapped into its old mesh. The vertex positions were determined 
by either of the two ways or by other prescribed methods. They were assumed to 
move with the fluid (Lagrangian mesh) 


*" +l = 


y?j +l = yfj 

(3.33.a) 


or they were not move with the fluid (Eulerian mesh) 

vtr = v® - (3 - 33 - 6) 

In Lagrangian technique, phase C is not needed. So in phase C, only the 
Eulerian calculation is discussed. From equation (3.33), the volume VT +1 of phase 
C was calculated by equations (3.1) ~ (3.4), the mass was calculated by equation 
(3.23) from continuity equation. 

M" +1 = M, B j +£)[/»]£ 6V <* 


(3.34) 
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where the 6V a is done by geometrically calculating the volume swept out by each 
regular cell face as ; t. moves from its Latrrmrgian position ( Jet.ertn.jneH by the .rf } . 
yfj) to its final position (determined by the . y?+ l )• the quantities in [p] are 
chosen from the i V a located. The relationships are shown in Figuie 3.5. 

Let SV a be the signed volume associated with the quadrilateral having swept b\ 
two opposing side s cell face o from phase B to phase n+1. This volume is evaluated 
by equations (3.1) ~ ( 3.4) and the algebraic sign of SV n can be decided. 

The density in n+1 phase was decided by 

yi + i _ v/ rl + 1 /V n_H (3.35) 

i J ij ~~ * 1 ij 


Similarly, the momentum equation was calculated by 

8 

Mtf'u ?? 1 = MfjU? + 

0-1 

8 


M n+ 1 v’'+ 

0 U ij 


1 ‘Vo 


3.36) 


where SV Q is the momentum cell similar to the 6V a for the regular cells. The a s in 
6V a is the a fac<‘ of the regular cell, a = 1, 2. 3. 4. The a s in <!>V 0 is the a face of 

the momentum ’ells, a = 1. 2. 3. 4. 5, 6, 7. S. 

The cell quantities [p]£, [/w]£, H« were determined by the donor cell scheme. 
Let the interested regular cell o be called cell 1 and the neighboring veil which is 
connected to the cell face o be called cell 2. The cell quantities were evaluated as 
an upwind-weigrted average of the quantity of cell face 1 and cell face _. First, the 
upwind cell or Conor cell should be decided by the sign of 8V : if 6V > 0. the cell 2 
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is the upwind or donor cell which cell 1 is the downwind or acceptor cell, and vice 
verse. 


Therefore, those quantities was defined 


nB _'Q 2 n sv a >o 

U ~ 1 Q? if 6V a < 0 

n B_ ( Q ? if 

_ 1 Q? if 6V a <0 


(3.37) 


where Q = p, pu.pt>, d is referred to the donor cell quantity and a is the acceptor 
cell quantity. In p]£, [pu]®, [pi>]£, a partial donor cell scheme was described. 




1 + «0 + doC) + ^Qa( 1 


0(1 


do C) 


;3.38) 


where o 0 and do are adjustable coefficients. (0 < a 0 + d 0 < 1) and 

mi 

V, + v 2 


(3.39) 


where Vj, Vi is t ie area of cell 1, cell 2. 

In equation (3.38), when o 0 = 0 and do = 0, this scheme represents the cen- 
tered differencina; approximation for the convection term, which is unconditionally 
unstable. If = : 1 and do = 0 , this scheme is the donor cell or upwind differencing 
scheme. It is stable, but too diffusive for most calculations. If o 0 = 0 and do = 
1, the scheme is an interpolated donor cell scheme. A weighted factor which is an 
average of the centered and the donor cell differencing is added, which makes this 
scheme more stable. This interpolation scheme is also less diffusive than the donor 


cell differencing by a factor C. 










Figure 3. 4 Regular Cell used for Differencing Viscous Stress 
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CHAPTER IV IMPLEMENTATION OF TURBULENCE 

MODELS IN THE ALE-ICE SCHEME 


The implementations of the two turbulence models described in Chapter 2 as 
well as the appropriate boundary conditions will be described in this Chapter. 


A. Subgrid Scale Eddy Viscosity 

The eddy viscosity fi t is a cell center variable and the velocity gradients in 
equation (2.8) ere evaluated for each of the two triangles into which the cell is 
divided (Figure 3.1) using the triangle rule: 


d<(> 

OX 


1/23 <jfrl3 ~ 2/13 023 
*13 J/23 - *23 1/13 


d<i> 

Ty )TR 


*13 </>23 ~ *23 </>13 
*13 2/23 — *23 2/13 


(4.1) 


where <p = u,v, ir w and 


fiab = 4>a - <t>b , 


X ab — X a i 


Vab — Ua Vb 


and 

a, b — 1 , 2 , 3, 4 

The derivatives of [d(j>/ dx) b l and (d(p/ dy) b L are obtained by the pei mutation 
(TR, 1, 2, 3) — ► (BL, 1, 3, 4). From those derivatives the quantities Dj } of equation 
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(2.9) was computed by 

Dfi — max[( Di-, )rp- {D tl )7*r ] 

Those results were substituted to equation (2.8), and the turbulent viscosity 
p t is obtained. 


B. fe - e Turbulence Model 

The k - e model requires solving two extra transport equations. In phase A, 
the Lagrangian formulation excluding the convection terms was implemented 


M A k A - M n k n _ A,, 


a 


fik 


G n - c D p n e r 


M A e A - M n e r 

At 


A f . fi .de M , 

= ^ iy 7, ) ai +{y ^ ) dy ) ° 


Of= 1 


(c u e n G n - c 2e p n e n )/k n 


(4.3) 


In phase C, the Donor cell scheme was used to calculated the convection terms 


C. Wall Functions 

The accurate solution of the boundary-layer equations for turbulent Hows using 
models which evaluate the turbulent viscosity at all points within the How lequnes 
that grid points be located within the viscous sublayer, y+ < 4.0 for incompressible 
How, and perhaps y + < 1.0 or 2.0 for flows in which a solution to the energy 
equation is also being obtained, where y+ = yu/v. The use of equal grid spacing 
for the transverse coordinate would require several thousand grid points across 
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the boundary layer for a typical calculation at moderate Reynolds numbers. This 
at least provides motivation for considering wavs to reduce the number of grid 
pants required to span the boundary layer, the techniques which have been used 
successfully fall into three categories, use of wall functions, unequal grid spacing, 
and coordinate transformations. The wall function will be disscussed in this section. 
For many turbulent wall boundary layers the inner portion of the flow appears to 
have a “universal’ character captured by the logarithmic “law of the wall”. The 
inner region is a zone in which convective transport is relatively unimportment. 
The law of the wall can be roughly thought of as a solution to the boundary 
layer momentum equation using Prandtl’s mixing-length turbulence model when 
convective and pressure gradient terms are unimportment. In this approach, the 
law of the wall is usually assumed to be valid in the range 30 < y + < 200 and the 
first computational point away from wall is located in this interval. The law of the 
wall equation was given by the logarithm form 

— = -ln ( — ) + £ (4 - 4 > 

U * K v 

where u is the tangential component of the fluid velocity at a perpendicular distance 
y from the wall, v is the molecular kinematic viscosity, k is von karman constantfin 
this study, the k is set to 0.4 ), B is a constant that depends on wall roughness (we 
set B = 5.5 (smooth-wall)), and u. is the shear speed, which is related to the wall 

shear stress t w and the fluid density p by r w = pu m . 

In this equation, the u. of equation (4.4) appears in both sides of the equation, 
which one would ordinarily have to solve iteratively. To avoid the inconvenience 



this would entail we use instead an approximation obtained by replacing u* in the 
argument of the logarithm by 1/7-law value u*. The 1/7-law’s equation (Hinze 
1975) is given by 

ii=8.3(^) 1/7 (4-5) 

It * V 

multiplied yu,/( 3.3i/) to both sides of equation (4.5) and took power 7/8, then the 
equation change to below 

^ = 0.15(^) 7/8 (4-6) 

V v 

substitute (4.6) into the right side of equation (4.4), the new equation was obtained 

— = 0.75 + 2.19 ln(— ) (4-7) 

u , v 

From this equation, it can be explicitly solved for u,. The wall stress is then 
given by t w = fu*. If yu/v (a Reynold’s number based on y ) is too small, the 
interested point lies in the laminar sublayer rather than the law-of-wall region. In 
this case, the laminar sublayer equation was given 

Ji = (^)i/2 (4.8) 

U» V 

The transition values between equation (4.7) and equation (4.8) is made at the 
point where they predicted the same u*, which is yu/v = 130.3. 

These equa ion are implemented numerically in the following way. Consider a 
typical cell adjacent to the wall, as shown in Figure 4.1. The shear stress t w for 


this cell is evaluated by setting 
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y = l(h L + h R ) ( 4 ‘ 9) 

and i. = //./,/p, where p is the density in the cell in question, and p L is evaluated at 
the temperature of the cell. Equation (4.7) and (4.8) incorporate two simplifying 
assumptions: (1) the normal velocity at points L and E is negligible, so that the 
tangential component may be replaced by the magnitude of the velocity, and (2) the 
lines connecting points L and R to their neighbors on the wall are approximately 
perpendicular to the wall. These values for u, y, and determine t*. and hence r*, 

as already described. 

With t w thuf determined, the product r w AAt gives the associated change m 
fluid momentum occurring on a time step (A is the wall area of the cell in question). 
Half of this change is apportioned to vertex L and the other half to vertex R. These 
changes are effect. *d by multiplying all velocity components at each of these vertices 
by a single factor for that vertex, so that the directions of the velocity vectors do not 
change. (This procedure again relies on assumption (T) of the previous paragraph.) 
These factors are given by 


Fl = 1-0 


r w AAt 


2M L (u 2 L + e|+iei)’ 


Fr = 1 - - 


,,AAt 


2 Mr(u 2 r + v\ + w 2 r )* 


(4.10) 


where M L and Ur are the vertex masses. If either of these factors comes out less 


than 0.3, It is replaced by 0.3 to prevent the velocities from changing by too large a 
factor on a single time step. It should be noted that these changes to the velocities 
at points L and R are only those due to the particular cell in question. A similar 
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change to the velocity at L will result from the wall cell to the left, and to the 
velocity at R from the wall cell on the right. 

The wall fur ction set for k is 


k = 



(4.11) 


and for e is 



ny 

and the turbulent viscosity /x* was calculated by equation (2.11). 


(4.12) 


D. Initial Conditions 

At the inlet of the calculation domain, all dependent variables have to be 
specified. The k profile may be estimated from the measurements. In general, k is 
specified as a percentage of the inlet mean square velocity. 

k in = 0.003 x u] n (4.13) 

The inlet profile for e has to be assumed, since no measurements are available. 
It is specified as 

f,„ = c^ 2 /(0.03D/2) (4.14) 

where D is characteristic length at the inlet and c ^ = 0.09. 

From this two initial values, the turbulent viscosity fit can be calculated by the 


equation (2.11). 




CHAPTER V RESULTS AND DISCUSSION 

For the purpose of code validation and testing, several benchmark problems 
with available analytical and experimental data have been chosen. The first case is 
the Burger s equation with complex geometries and boundaries. Next, the laminar 
and turbulent pipe flows are calculated. Finally, a confined coaxial jet flow relevant 
to combustor des gn is calculated. 


A. Solution Of Biir^er’s Equa tion With 
Complex Qpomotries And B oundaries 

The purpose of this problem is to test the capability of the ALE-ICE scheme for 
handling complex geometries and boundaries. The Burger’s equation with mixed 
Neumann and Dirchlet boundaries are chosen to test the ALE-ICE scheme. 

The Burger s equations for the unsteady incompressible flow may be written 


as 


where 


du du du d 2 u d u r 

-p U “ — "P — t/( r\2 T 09 / W Jx 

dt dx dy d 2 x 


d 2 y 

d 2 v 


dv dv dv d 2 v , 

di + u di + V dy + 3 2 y 


)+/, 


/, = + 3*V - Ivy 

/v = rT7^+2xy-^) + 3xV-2- 


(5.1) 
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with the Neumann boundary conditions and the Dirchlet boundary conditions spec 
ified according to Figure 5.1. The Neumann boundary conditions ax* soecified: 


du 

dx 


2 xy 


du 

dy 


= x 


2 


dv 



2 



The solution will be compared with the exact solution 


(5.2) 


u 


1 

r+7 


+ x 2 y 


v 


1 +t 


+ xy 


(5.3) 


First, grid independence studies were carried out by using three grids . 4 x 
10, 8 x 20 and 16 x 40. The case considered was v = 1 and times step was set rather 
small to reduce any instability. The total r.m.s. errors, which was defined as sum 
of the pointwise r.m.s. error 

[W = {£<„,. -s.)V£»7) ,/2 1 < 5 - 4 > 

were plotted vs time on Figure 5.2. Both 8 x 20 and 16 x 40 reaches asymp- 
totic values of 7%, and since only the long- time behaviors are compared the rest 
of the results are shown using the 8 x 20 grid. Effect of viscosity on the calcu- 
lated results is shown on Figure 5.3. As expected, as Reynolds number increases 
( v decreases), convection term dominates and due to the upwind scheme for the 
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convection term approximation, the total r.m.s. errors increases. Typical results m 
terms of pointwise errors are shown in Figure 5.4 and Table 5.1. For areas closed 
to non-orthogonal grids and Neumann's B.C.'s. there appear to have less accurate 
results (Lee et al. 1989). This points to the further development in terms of flux 
calculations in the ALE-ICE method for irregular grids and Neumann Boundary 
condition implementations. 


B. Laminar P ipe Flow 

The second testing case is a laminar pipe flow. The pseudocompressibility 
method describee in Chapter 3 was used for the incompressible flow case. This 
method has also been used in Rogers et al. ( 1987), Liang et al. (1985) and Nichols 
et al. (1980) for incompressible flows. 

The analytical solutions of laminar pipe flows were given in Bird et al. (1960). 
The calculated results will be compared with the analytical solution. The velocity 
distribution in the fully developed region was given by 

C f,=2C/s|l-(|) 2 ] < 5 ' 5 > 

where U. is the .ongitudinal velocity. U t is the bulk velocity and R is the radius of 
the pipe. 

The Reynolds number (Re = DU b /u) of 1,000 was carried out for laminar flows 
in which D is the pipe diameter and v is the viscosity. In Figure D.o, 0.6 different 
3 values (pc 2 ) were tried to find the optimal performance. For 3 around 5 the 
solutions obtained were optimal. Next, grid independence studies were earned out 
by using 11 x 81, 21 x 81 and 31 x 81 grids. In Figure 5.7, 5.8 fully developed 
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and developing profiles show that 21 x 81 and 31 x 81 grids converged to the same 
solutions end compare well w.th the analytical solution velocity within about .8 % 
errors. The calculation was done by assuming the plug profile at the pipe inlet and 
the calculation domain was 60 in pipe diameter. Fully developed (Figure 5.9) and 
developing (Figure 5.10) data compared to analytical solution velocity is only .75 % 
errors in 21 X 81 case. In this section, each case runs about 17,000 cycles and each 
cycle uses .7 second on CRAY-XMP. This simulation will give some suggestion to 
the ALE-ICE scheme's future works and will help the next test case for turbulent 

flows. 


C. Turbulent Pipe Flow 

Most turbulence models have been tested and/or developed using fully devel 
oped pipe flow conditions, and this testing will be made m this study. As recom- 
mended by Mari inuzzi and Pollard (1989), the calculated results are compared to 
one another, and, for the sake of placing the calculations in perspective, they are 
compared to various data sets that appear in the literature. A single data set has 
not been used because, in the authors’ opinion, a reliable, well- documented data 
set does not exist due to back of providing detailed inlet conditions. As a results, 
the data sets o’ Richman and Azad (1973), Lawn (1971) and Nikuradse (1932) 
will be compared with the calculation results. From Martinuzzi and Pollard ( 1989) 
recommendation, the Re =10,000 and Re =100,000 will be used in this study. 

Calculation grids 21 x 81, 31 x 81 for Re =10,000 and 51 x 81. 61 x 81 and 71 x 
81 for Re =100 000 were used by grid independence study. The plug flow condition 
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was used at the inlet. The (3 value for 21 x 81 and 31 x 81 was 500 and for 51 x 
81, 61 x 81 and 71 x 81 grids was 5 000. The calculation domain was 80 in the 
pipe diameter. Empirical relationships are used to assign entrance values to k and 
e; that is, k = 0 )03 U% for Re =10.000, k = 0.03 U* for Re = 100,000 and 6 = 

(C„Jfct/0.03/J*). 

In Figure 5.1 1, the fully developing axial velocity profiles normalized by the in- 
let velocity, at Re =10,000, was presented. Figure 5.12, the fully developed velocity 
profiles are compared with those data of Nikuradse. The fully developed profiles 
at Re =10,000 for u + = u/u* vs y + is shown in Figure 5.13, the slopes are gener- 
ally well reproduced, but the absolute values of the calculations tend to overpredict 
the data of Nikuradse (1932). The developing velocity profiles for Re =100,000 
are shown in Figure 5.14 with the data of Richman and Azad (1973). The differ- 
ences are negligible. The fully developed axial velocity profiles at Figure 5.15, at 
Re =100,000. aie compared with those data of Richman and Azad (1973). The 
agreement between experimental data and calculations is poorer close to the wall 
boundary. This c an be further seen from Figure 5.16 where the mean velocity profile 
is plotted using he wall coordinate. Using the 1/7-th law wall function forced the 
first grid calculation results to match the logarithm wall. Due to the poor turbu- 
lence model capability for strong shear flows, the predictions deviate significantly 

away from the wall. 

The distribution of the turbulence kinetic energy in the fully developed region 
for Reynolds number of 10,000 is overpredicted (Figure 5.17) and is underpredicted 
(Figure 5.18) for- Re =100,000 as compared to Lawn’s experimented data. 
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It should be noted that the 1 /7-th law wall function used in there calculations 
is otdv valid for fully developed turbulent pipe flows. This assumption is incapable 
of modeling developing regions and separated regions as will be seen in the next 
test case. The use of the 1/7-th law wall function is probably responsible for the 
poor performance of the k - e model toward the wall regions as implemented in this 

study. 

D. The Coaxial Jet Flow 

The Space Shuttle Main Engine (SSME) high pressure fuel prebumer consists 
of 264 coaxial injector elements injecting gaseous hydrogen and liquid oxygen into 
a large cylindrical chamber. The ARICC code was originally developed to simulate 
these types of flowfield. The outer wall is a free-slip boundary that corresponds 
to the imaginary symmetry surfaces of a streamtube whose diameter is determined 
based on the average cross-sectional area available to each injector element (Liang 
et al. 1986a). This ability is deemed necessary for the simulation of the confined 
coaxial jet flowt. Owen’s turbulent recirculating jet flows (1976) will be tested in 
this study. 

Owen’s turbulent recirculating jet flows (1976) is a complex one consisting of 
large recirculating zones. It is a confined coaxial jet flow with high velocity ratio 
between annular jet and central jet, and a central toroidal recirculation zone is 
formed due to the imbalance of the mass flow between the central and annular jets. 
The sudden expansion of annular jet results in a severe adverse pressure gradient. 
When this pressure gradient is too strong for the low momentum central jet, a 
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central toroidal recirculation zone results in addition to the corner recirculation 
zone. 

In an effort o provide such data and also to understand the nature of the 
flow in a coaxial jet flow configuration of interest to combustor designers, Owen 
carried out detaile d turbulence measurements using LDA. The Owen s experimental 
configuration (reproduced in Figure 5.19) consist of a 2.5 in (6.35 cm) central jet 
surrounded by a 1.5 in (8.89 cm) annular jet. The jets discharge into a o m (12.7 
cm) diameter chamber 4S in (121.9 cm) long. The outer and inner peak velocities 
were 96.0 (29.26) and 8.0 (2.44) fps (m/sec) corresponding to Reynolds numbers 
based on respective jet diameters of 1.5 and 0.08 x 10 5 . Inflow profiles were not 
reported. 

Sturgess, et al. (1983), Syed and Sturgess (1980), Novick. et al. (1979) pre- 
dicted the confined coaxial jet flow corresponding to the experimental conditions of 
Owen’s (1976) w:th k-e turbulence model. None of these predictions produced the 
shape, size, and location of the important central toroidal recirculation zone cor- 
rectly (Nallasamv 1987). The length of the central recirculating zone was severely 
underpredicted by about 40 percent. As concluded by Nallasamy (1987), the corner 
and central recirculation zones are very sensitive to the central jet exit geometry. 

Computations for this case were made using ARICC code with 21 x 81 and 31 
x 81 unstructured grids, shown in Figure 5.20. The calculation domain was 12 times 
of the diameter of the chamber and extended 4 inches for inlet flow. The inner jet is 
2.5 inches, the cuter jet is 3.5 inches and the chamber is 5.0 inches. The Reynolds 
number for inner jet is 8000 and for outer jet is 15,000; Re = DUb/v , where D 
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is the jet correspondence diameter and U b is the bulk velocity for correspondence 
jets, v is the kinematic viscosity. The inlet values for k is 0.003 xT\ 2 and for e is 
/0.03/.R. C imputations were made with the subgrid and the k - e model. 

Figure 5.21 shows the comparisons of the centerline mean axial velocity profiles 
and the measured values. The k - e model prediction is good up to 2.0 D and the 
subgrid model only within 0.5 D. Beyond those the predicted value are smaller than 
the measured ore. The mean velocity profiles predicted by the k - e model and 
subgrid model compared to the measured values are illustrated in Figure 5.22, 5.23. 
The predicted mean velocity profiles of k-e model are m good agreement with the 
data, but subgrid model are not. Between x/D=0.25 and x/D=1.8, some differences 
between the data and the predictions are seen. However, the overall predictions of 
k - e model are in good agreement with data in the recirculation zone. The k - e 
model predictions for r.m.s. velocity fluctuation and the data are shown in Figure 
5.24. The prediction is very close to the data. Figure 5.25 shows the streamline 
patterns, subgrid model center recirculating zone is larger than that of k-e model s. 
The k-e model shows a better representation of the experimental data. For corner 
recirculating zones, k — t model also shows good predictions 

Figure 5.26 shows the transient behavior of ARICC code. As seen from those 
profiles, they are not real transient results. Equation 3.22 indicated that the con- 
tinuity equation can only be satisfied steady state was reached From this point of 
view, the transient results are not time-accurate and only represent intermediate 
transient to the final steady state results. 
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Figure 5. 1 Complex Geometries Grid used in Burger’s Equations 
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Figure 5. 2 Total R.M.S. Error for Grid Independent Study 
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Figure 5. 4 8 x 20 Grid Cell 
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Table 5. 1 Values for Figure 5.4 



error 

rms 

u 

V 

4.9E-5 

9.18E-5 

1.12E-4 

2.73E-5 

1.85E-4 

9.06E-5 

6.7E-5 

9.6E-4 

4.6E-3 

8.9E-3 

6.4E-3 

1 .06E-2 

3.4E-3 

9.9E-3 

1.24E-3 

4.24E-3 

7.61E-4 

2.28E-3 

6.2E-4 

2.18E-3 
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Figure 5. 5 Fully Developed Velocity Profiles with Different /I Number 
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Figure 5. 6 Development of Centerline Velocity for Different 0 Values, 0-5, 50, 500 
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Figure 5. 8 Axial Velocity Profiles for Normalized Data 
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Figure 5. 9 Fully Developed Velocity Profiles 
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Figure 5.10 Axial Velocity Profiles 
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Figure 5.11 Development of Centerlines Velocity Profiles for Re =10,000 
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Figure 5.12 Axial Velocity Profile for Re =10,000 
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Figure 5.13 vs at y/D=80 for Rc —10,000 
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Figure 5.14 Development of Centerline Velocity Profiles for Re =100,000 
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Figure 5.15 Variation of Axial Velocity with y/D=40, 50, 70 for Re -100,000 
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Figure 5.16 u + vs y + at y/D=80 for Re -100,000 
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Figure 5.17 Turbulence Kinetic Energy vs x/D at y/D-80 for Re -10,000 
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Figure 5.18 Turbulence Kinetic Energy vs x/D at y/D-80 for Re -100,000 
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Figure 5.20 Unstructure Grid for Confined Coaxial Jet Flow 
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e 5.22 Axial Velocity at y/D-0.25, 0.6, 1-0 
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Figure 5.23 Axial Velocity at y/D-1.4, 1.8 
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Figure 5.24 Axial Turbulence Intensity on the Centerline 
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Figure 5.26 Transient Delta. 
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CHAPTER VI CONCLUSION 


In this study, the ALE-ICE scheme has been applied to four different case, the 

solved on a complex geometries with different boundary 

solved in other cases. The pseudo- 


Burger’s equations were 

conditions. Fully Navier-Stokes equations were 

used to solve incompressible fluid flow problems. Both 


tested and compared with the analytical so- 


compressibility me', hod was 
laminar and turbu.ent pipe flows were 
lution (Bird et al. 1960) and experimental data and other predictions (Martinuzri 

et al. 1989). Last y, a confined coaxial jets flow with larger momentum difference 
between inner je, and annular jet were calculated. Owen’s experimental data was 

chosen for comparison. 

Results obtained from solutions of the Burger’s equation by unstructured grid 
showed that the lotal r.m.s. errors, defined by equation (5.4), reached asymptotic 
values of 7%. Effects of viscosity on the calculated results were studied. As ex 
pected, as Reynolds numbers increase („ decrease), convection term dominates and 
the total r.m.s. errors increase. For areas closed to non-orthogonal grids and Neu- 
mann B.C., there appear to have less accurate results. This points 
development in terms of flux calculations in the ALE-ICE method for irregular 
grids and Neumann Boundary condition implementations. Results from laminar 
pipe flow showed good agreement with the analytical solution (Bird et al. I960). 
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From turbulent pipe flows results, the calculations using the i - « -ode. showed 
good predictions except the near wall regton. the discrepancies were traced to the 
usage of the l/7th law wall function in which the fully developed assumptions 
used. For the confined coaxial jet case, the subgrid model can not matched Owen's 
experimental data 1976). on the other hand, the h - c model gave the closer pre 
dictions compared to the experimental data. For these tests, the adiabatic sound 
speed factor, 0 = for the pseudocompressible continuity equation can not be 

obtained easily. The number of d suggested by Nichols et al. ,1980). Anderson et 
al. (1984), Soh (1987) and Rogers (1987) could not be used in ARICC 
different numerical method. Liang (1985) suggested that 0 could be any suitably 
large number to convert the small volume changes into appropriately large pressure 
changes. However, the 0's being used in laminar pipe flows were found to be dif- 
ferent from those of the turbulent pipe flows. Also, the 0 value used in the subgrid 
model calculations was different with those used in the t ~ « model calculations. 

It is concluded that there is no specific 3 that can be used in any cases. Trial 
and error has to be executed for each different cases. Fhom this point of view, the 
pseudocompressibility method is more difficult for general fluid flow simulations. 

From the results of confined coaxial jet calculations (Figure 5.27), the transient 
behaviors of the pseudocompressibility method in the ARICC code were not time- 
accurate. Those profiles look like a temporary acoustic waves pushing through the 
calculation domain. It is similar to the compressible flows calculations. Therefore, 
the ARICC cone is not capable of time-accurate predictions of the liquid jet in 
two-phase confined coaxial jet simulations. In this regard, a subiteration scheme 
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can be incorporated in pseudotinre to satisfy the continuity equation at each time 

step as done recently bv Rogers snd Kwak (19^0). 

From the cases tested in this study, it is aiso concluded that the ARICC code 

using the pointw.se SOR solver is not efficient for fine grid calculat.ons. A more ef- 
ficient solver such a. the conjugate gradient method should be implemented for the 
future usage. As for the turbulence models, the more elaborated k - s model gener- 
al gave more satisfactory predictions in the cases tested compared to the subgrid 
model. The .mplen.entation of the l/7th law wall function for the near wall velocity 
profiles was found to Ire .naccurate for developing flows. Further improvements on 
the near wall treatments of turbulence models are definitely needed. 
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